Thermodynamic approach to dense granular matter: a numerical realization of a 

decisive experiment. 
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Years ago Edwards proposed a thermodynamic description of dense granular matter, in which the grains (the 'atoms' 
of the system) interact with inelastic forces. The approach is intriguing but is not justified from first principles, and 
hence, in the absence of conclusive tests of its validity, it has not been widely accepted. We perform a numerical 
experiment with a realistic granular matter model specially conceived to be reproducible in the laboratory. The results 
strongly support the thermodynamic picture. 



The similarities between a driven granular system and 
a system of molecules in a fluid have prompted the use of 
statistical ideas to describe granular matter |l|, such as 
for example kinetic theories of gases to study rapid flows 
of low-density systems of inelastic particles |^. In this 
context, it is customary to associate the mean kinetic en- 
ergy of the particles with a 'granular temperature', hav- 
ing no real thermodynamic meaning. 

For the opposite limit of dense, slowly flowing gran- 
ular matter, a more radical train of ideas was initiated 
more than a decade ago by Edwards and coworkers: the 
proposal of a statistical ensemble y-§|, and through it, 
thermodynamic notions such as entropy and tempera- 
ture (the latter unrelated to the 'granular temperature' 
above). Although the idea of a thermodynamic descrip- 
tion of granular matter was recognized as attractive, it 
was not universally accepted because there is no known 
first principle justification of Edwards' statistical ensem- 
ble, such as there is for ordinary statistical mechanics of 
liquids or gases (Liouville's theorem). 

Recently, analytical developments and numerical work 
on schematic models originally devised for glasses have 
given a new perspective for the dynamics of dense gran- 
ular matter |7|-p^. The theoretical support provided by 
this framework for a statistical approach has spurred a 
renewed interest in Edwards' thermodynamics. 

On the experimental side, years ago Nowak et al. ||lj] 
studied in detail the density fluctuations in a vibrated 
granular material, and proposed that these fluctuations 
should reflect an underlying thermodynamics, much as 
they do in an ordinary thermal system. However, as far 
as the actual verification of Edwards' approach, the ev- 
idence they found was if anything rather negative: we 
shall discuss below some possible reasons for this. On 
the other hand, Edwards' approach has never been tested 
with simulations of realistic models of granular materi- 
als, characterized by the fact that energy is supplied by 
external driving (via shear or tapping) and dissipated by 
inelastic collisions and slippage between the grains. 

In this work we describe a numerical experiment spe- 
cially conceived to be reproducible in the laboratory, us- 
ing a realistic model of sheared granular matter. Con- 



sidering particles of different sizes inside the medium, we 
compute a dynamical temperature from an Einstein rela- 
tion relating random diffusion and mobility (transverse to 
the shear direction) at long time scales. If there is an un- 
derlying thermodynamics, it will impose that the dynam- 
ical temperature be independent of the tracer's shape: a 
very strong condition which we show to hold. We then 
perform an explicit computation to show that the tem- 
perature arising from Edwards' thermodynamic descrip- 
tion and the dynamic one measured from Einstein's re- 
lation coincide. This last step cannot be performed in 
the laboratory, so the numerical simulation provides the 
missing link between thermodynamic ideas and diffusion- 
mobility checks. 

Consider a 'tracer' body of arbitrary shape immersed 
in a liquid in equilibrium at temperature T. As a con- 
sequence of the irregular bombardment by the particles 
of the surrounding liquid, the tracer performs a diffusive, 
fiuctuating 'Brownian' motion. The motion is unbiased, 
and for large times the average square of the displace- 
ment goes as (|x(i) — x(O)p) = 2Dt, where D is the 
diffusivity. On the other hand, if we pull gently from 
the tracer with a constant force /, the liquid responds 
with a viscous, dissipative force. The averaged displace- 
ment after a large time is {[x{t) — x(0)]) = fxt, where 
X is the mobility. Clearly, the same liquid molecules are 
responsible both for the fluctuations and for the dissipa- 
tion. Although both D and x strongly depend on the 
shape and size of the tracer, they turn out to be always 
related by the Einstein relation D/x — T, where T is the 
temperature of the liquid. 

Conversely, if in a fluid of unknown properties we flnd 
that several tracers having different diffusivities and mo- 
bilities yield the same ratio D/x = T, we may take this 
as a strong evidence for thermalisation at temperature 
T. Indeed, recent analytic schemes for out of equilibrium 
glassy dynamics have suggested the existence of such a 
temperature T governing the slow components of fluctu- 
ations and responses of all observables ^j^ ■ 

In order to test the existence of this temperature for 
dense granular matter, we perform, with a realistic nu- 
merical model, a diffusion-mobility experiment in con- 



ditions that can be reproduced in the laboratory. In 
the model, deformable spherical grains interact with 
one another via non-linear elastic Hertz normal forces, 
and non-linear elastic and path-dependent Mindlin trans- 
verse forces pj]. The normal force, F„, has the typical 
3/2 power law dependence on the overlap between two 
spheres in contact, while the transverse force. Ft, de- 
pends on both the shear and normal displacements be- 
tween the spheres. For two spherical grains with radii 
Ri and R2: Fn = I k^R^/^w^/^, AFt = kt{Rw)'^/^As. 
Here R = 2RiR2/{Ri + R2), the normal overlap is 
w = (1/2) [(i?i -I- R2) — \xi — X2\] > 0, and xi, X2 are 
the positions of the grain centers. The normal force acts 
only in compression, _F„ = when w < Q. The vari- 
able s is defined such that the relative shear displace- 
ment between the two grain centers is 2s. The prefactors 
kn = 4G/(1 — v) and kt — 8G/{2 — v) are defined in 
terms of the shear modulus G and the Poisson's ratio v 
of the material from which the grains are made (typically 
G = 29 GPa and v = 0.2, for spherical glass beads). We 
assume a distribution of grain radii in which Ri = 0.105 
mm for half the grains and R2 = 0.075 mm for the other 
half. The observables are measured in reduced units: 
length in units of R, force in units of GR^ , time in units 
of ^JpR?/G, where p is the density of the particles. Inter- 
nal dissipation is included in two ways: (1) via a viscous 
damping term proportional to the relative normal and 
tangential velocities, and (2) via sliding friction: when 
Ft exceeds the Coulomb threshold, /iF„, the grains slide 
and Ft — fJ-Fn, where fi is the friction coefficient between 
the spheres (typically /i — 0.3). 

We perforin molecular dynamics (MD) simulations for 
a binary system of 200 large and small spheres in a pe- 
riodic 3D cell. Our calculation begins with a numerical 
protocol designed to mimic the experimental procedure 
used to prepare dense packed granular materials. The 
simulations start with a gas of spherical particles located 
at random positions in a periodically repeated cubic cell 
of side L. The system is then compressed and extended 
slowly until a stationary situation with a specified value 
of the pressure and volume fraction — above the random 
close packing fraction — is achieved. We then apply a 
gentle shear in the y-direction at constant volume by 
moving the periodic images at the top and bottom of the 
cell with velocities jL/2, where 7 is the shear rate (Lees- 
Edwards boundary conditions) . In this simple shear flow 
the gradient of the velocity along z is uniform (no shear 
bands). We apply a shear rate of 7 = 10^**, which allows 
the system to be in the quasi-static regime. We focus 
our study on the slow shear rate, quasi-static limit where 
the system is always close to jamming, but moves just 
barely enough to avoid stick-slip motion fia]. Here the 
external pressure is large enough ('^ 10 MPa in a typical 
simulation) and the average coordination number is high 
enough (typically '^ 7) that deformation and elasticity 
of the particles play the dominant role, as opposed to 



the collision dominated rapid flow regime described by 
kinetic theories. We checked that shear induced segrega- 
tion is absent at the times scales of our simulations. The 
contacts between particles are enduring and the internal 
stresses in the system are transmitted via a network of 
'force chains' M,0|. independent of the shear rate. 

After a transient of the order of the inverse shear rate, 
we start measuring the spontaneous {\x{t) — x{0)\'^) and 
force-induced displacements {[x{t) — x{0)])/f along the 
x-direction for the two types of particles with different 
sizes. The results are shown in Fig. fl. We notice that the 
diffusivities and the mobilities are different for the two 
type of particles, since they have different sizes. However, 
when we draw the parametric plot of {\x{t) — x(0)|^) ver- 
sus {[x{t) - a;(0)])// (Fig. ||) we find parallel straight 
lines for large time scales, implying an extended Einstein 
relation: 
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valid for both particles with the same Tdyn for widely sep- 
arated time scales t ^ tw This suggests that Tdyn can 
be considered to be the temperature of the slow modes. 

For small time scales (fast rearrangements, near the 
origin in Fig. m the fluctuation-dissipation plot shows 
that there is not a well defined Einstein-relation tem- 
perature. We expect this result since the fast motion of 
the grains depends on the microscopic interactions dom- 
inated by inelastic collisions |0]. We also calculate the 
'granular kinetic temperature' defined in terms of the ve- 
locity fluctuations in the a;-direction. Unlike Tdyn, we 
find that this temperature is different for the two types 
of particles, and their values are two orders of magnitude 
smaller than Tdyn- We expect this result since the kinetic 
granular temperature is dominated by the fast (high fre- 
quency) modes, which are not thermalised to Tdyn- 

We also repeat the numerical experiment for a sys- 
tem of Hertz spheres without transverse forces (n — 0: 
experimental realizations of elastic spheres with viscous 
forces but without sliding friction are foams and com- 
pressed emulsions |l8|Jl9|1 ) and find that Tdyn is well de- 
fined at long time scales for this case as well (see Fig. |^). 
Thus, our results suggest that the validity of a structural 
temperature for long-scale displacements (larger than a 
fraction of the particle size) holds in the presence of vis- 
cous forces between grains or even of a sliding threshold 
(Coulomb's law). 

The existence of a single temperature is an instance 
of the zero-th law of thermodynamics, for which we find 
positive evidence here, at variance with the experimen- 
tal result in Nowak et al ]13] . Three possible reasons for 
the apparent violation of the 'zero-th law' in their exper- 
iments are i) the effect of an unknown height-dependent 
pressure (as pointed out by authors), ii) a rather high 
tapping amplitude (Edwards ensemble is in principle only 
valid in the high compaction limit pX|]) and iii) the fact 



that the density fluctuations considered are integrated 
over aU frequencies, thus including also 'fast' relaxations. 
Next, we treat the question whether it is possible to 
relate the dynamical temperature obtained above to the 
thermodynamic construction of slowly driven out of equi- 
librium systems proposed by Edwards. Whereas in the 
Gibbs construction for equilibrium statistical mechanics 
one assumes that the physical quantities are obtained as 
average over all possible configurations, Edwards ensem- 
ble consists of only the blocked configurations (static or 
jammed) at the appropriate energy and volume. The 
strong 'ergodic' hypothesis is that all blocked config- 
urations of given volume and energy can be taken to 
have equal statistical weights. This formulation leads 
to an entropy SEdwiE,V), and the corresponding tem- 
perature T^J^ = dSEdw/dE and compactivity ^^]„ = 

dSEdn^/dV. 

In order to calculate Tsdw and compare with the ob- 
tained Tdyn we need to count the number of blocked con- 
figurations at a given energy and volume. (For this cal- 
culation we concentrate in the case without tangential 
forces and sliding friction, in order to avoid path depen- 
dency which would lead to an ambiguity in the defini- 
tion of blocked configurations — see below). To do this 
in practice, we extend the 'auxiliary model' method fin] 
to the case of deformable grains. We define an auxiliary 
model composed of the true deformation energy {E, the 
Hertzian energy of deformation of the particles) and of 
a term that vanishes in (and only in) the blocked con- 
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(Here Fa is the total 



contact force exerted on particle a by its neighbours). 
We perform equilibrium MD with two auxiliary tempera- 
tures (T*,Taux), corresponding to the partition function 
X]exp[-i?/r* - Ebiock/Taux]- Anneafing Tai^x to zero 
selects the blocked configurations [Euock — 0), while T* 
fixes the energy E. We start by equilibrating the system 
at high temperatures {Taux and T* ~ oo) and anneal 
slowly the value Taux to zero and tune T* so as to reach 
the observed value of E during shear. At the end of 
the process {Taux — > 0), T*{E) — TEdw{E), since in this 
limit we are sampling the configurations with vanishing 
fraction of moving particles at a given E. 

Figure ^ shows various annealing protocols. We plot 
the value of the true deformation energy as a function of 
Taux during annealing using different values of T* as in- 
dicated in the figure. At the end of the annealing process 
the system stabilizes at a given deformation energy. Only 
when T* is equal to Tdyn obtained through the Einstein 
relation during shear, we find that the final E coincides, 
within the accuracy of the simulations, with the mean 
energy of the system under shear. For other values of T* 
the annealed system equilibrates at energies which are 
above of below the distribution of deformation energies 
during shear as seen in Fig. S. This proves the equality 
between Edwards' and the dynamic temperature. 



We conclude with some remarks: i) Since the blocked 
configurations are the same whatever the inter-grain dis- 
sipation coefficient, Edwards' ensemble (and hence its 
temperature) are insensitive to viscous dissipation. 

ii) On the contrary, tangential forces and sliding fric- 
tion block certain configurations, depending on how they 
are accessed. The ensemble of blocked configurations is 
then ill-defined. We have not tried to construct any en- 
semble, but content ourselves with the observation that 
Tdyn is also in this case independent of the particle size 
(Fig. |2|) — we suspect that thermodynamic concepts ap- 
ply, but that the relevant ensemble goes beyond Edwards' 
construction as it stands. 

Hi) We have tested the validity of the thermodynamics 
in an ideal homogeneous system with periodic boundary 
conditions by explicitly avoiding structural features of 
dense granular flows such as shear bands and segregation 
of the species. Even though it remains to be seen whether 
the thermodynamic picture can account for these inho- 
mogeneous effects, our ideal system may prove to be use- 
ful in deriving constitutive relations to be used in macro- 
scopic theories of slow granular flows. 

To summarize: we have performed numerically an ex- 
periment with dense granular systems specially conceived 
to be a 'dress rehearsal' for the real laboratory one. 
The independence of the Einstein-relation temperature 
on the tracer provides a strong test for the thermody- 
namic ideas. We have also showed that this temperature, 
obtained from measurable quantities, indeed matches Ed- 
wards' one. If, as now seems likely, this result is con- 
firmed in the laboratory, this will give an experimental 
foundation for the use of the powerful tools of statistical 
mechanics as the framework to study this kind of far- 
from-equilibrium dissipative dynamical system. 
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FIG. 2. Parametric plot of diffusion vs response function 
for small and large grains and for spheres interacting with 
tangential forces and without tangential forces (Coulomb fric- 
tionless). The fitting at long time scales shows the existence 
of a well defined temperature which is the same for small 
and large grains: Tdy„ = 2.8 x 10~^ for grains without trans- 
verse forces and Tdyn ~ 1.2 x 10"* for grains with Mindlin 
transverse forces and Coulomb friction. We calculate the re- 
sponse function for several small external fields and find the 
same temperature indicating that we are in the linear response 
regime. Plotted are results for a system without transverse 
forces using: / = 1.7 x 10~^ (small grains D, large grains o) 
and / = 2.6 x 10~^ (small grains ■, large grains •). For a 
system with tangential forces and Coulomb friction we show 
the case / = 6 x 10~^ (small grains A, large grains x). 
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FIG. 1. (A) Diffusion and (B) response function for the 
large and small particles in a sheared granular material mea- 
sured perpendicular to the shear plane as a function of time 
in MD steps. Both quantities depend linearly on time at the 
late stages of the evolution. The response function is mea- 
sured by applying a constant force / in the x-direction to 
each type of particle. Averages are taken over 30 temporal 
realizations. The obtained diffusivities and mobilities depend 
on the particles size as expected. 
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FIG. 3. Annealing procedure to calculate Tsdw at dif- 
ferent true energies. We plot the true deformation energy 
versus Taux together with the distribution of deformation 
energies obtained during shear (dashed curve, mean value 
(E) = 8.4 X 10""'). We equilibrate the system for 40 mil- 
lion iterations at A: (T* = 3.4 x 10~^,Taux = 3 x 10"*). 
We then anneal slowly both temperatures until B: 
(T* =3.4xl0"'',r„„a; =3x10"^°), where we split the trajec- 
tory in three paths in the {T* ,Taux) plane. Path 1: we anneal 
Taux -^ and T* — > 2.8 x 10~^ which corresponds to Tdyn 
obtained during shear (Fig. bf). Path 2: we anneal Taux -^ 
and T* ^ 3.4 x lO"*^. Path 3: we anneal Taux -^ but keep 
T* = 3.4 x 10"'* constant. When we set T* = Tdyn (Path 1), 
the final true energy value when Taux —> falls inside the dis- 
tribution of energies obtained during shear and it is very close 
to (E). This proves that Tdyn = Tsdw under the numerical 
accuracy of the simulations. For other values of T* 7^ Tdyn 
the final E falls out of the distribution obtained during shear 
(Paths 2 and 3). We also follow different trajectories (not 
shown in the figure) to T* -^ 2.8 x 10~^, Tau^, -^ and find 
the same results indicating that our procedure is independent 
of the annealing path. 



